import numpy as np
import pandas as pd
import math
import json
from matplotlib import pyplot as plt
from scipy.integrate import quad
from datetime import datetime, timedelta
# @params:
#     SAModel:需要添加该约束的模型
#     n_AC:空调数量
#     Ca:气体热容序列
#     Ra:气体热阻序列
#     Cm:气体热容序列
#     Rm:气体热阻序列
#     Ta0:室内气体初始温度序列
#     Tm0:室内固体初始温度序列
#     Pmax:最大制冷功率
#     Tmax:室内气体最高温度
#     Tmin:室内气体最低温度
#     To:室外温度序列
#     Cope:能效比序列
#     periods:时间段
def TempGenerator(periods):

    # 温度模拟（根据规律生成数据）
    # 清晨(00:00-06:00)温度缓慢上升
    morning_temps = np.linspace(22, 26, round(periods/4))

    # 上午(06:00-12:00)温度快速上升
    noon_temps = np.linspace(26, 32, round(periods/4))

    # 下午(12:00-18:00)温度维持高位
    afternoon_temps = np.linspace(32, 30, round(periods/4))

    # 晚上(18:00-00:00)温度下降
    evening_temps = np.linspace(30, 25, round(periods/4))

    # 合并所有时间段的温度数据
    temps = np.concatenate((morning_temps, noon_temps, afternoon_temps, evening_temps))

    return temps

def ACPowerGenerator(periods,n_AC,random_seed=40):
    ACParam = {}
    step = 24 / periods
    np.random.seed(random_seed)
    Ca = np.random.uniform(0.09, 0.11, size=n_AC)
    Cm = np.random.uniform(1.4, 1.6, size=n_AC)
    Ra = 0.2*np.random.uniform(18, 27, size=n_AC)
    Rm = np.random.uniform(1.8, 2.7, size=n_AC)
    Cope = np.random.uniform(2.8,3.2, size=n_AC)
    Ta0 = 26*np.ones(n_AC)
    Tm0 = 26*np.ones(n_AC)
    Ta = 26*np.ones((n_AC,periods))
    Tm = np.zeros((n_AC, periods))
    Pac = np.zeros((n_AC,periods))
    To = 1.2*TempGenerator(periods)

    for i in range(n_AC):
        # 计算解析递推形式的参数
        b = 1 / (Ca[i] * Ra[i]) + 1 / (Ca[i] * Rm[i]) + 1 / (Cm[i] * Rm[i])
        c = 1 / (Ca[i] * Cm[i] * Ra[i] * Rm[i])
        delta = b ** 2 - 4 * c
        if delta > 0:
            A = np.zeros((periods, 2, 2))
            B = np.zeros((periods, 2))
            r1 = (-b + math.sqrt(b ** 2 - 4 * c)) / 2
            r2 = (-b - math.sqrt(b ** 2 - 4 * c)) / 2

            e1 = math.exp(r1 * step)
            e2 = math.exp(r2 * step)
            for k, t in enumerate(np.arange(0, 24, step)):
                A[k, 0, 0] = (r1 * e1 - r2 * e2) / math.sqrt(delta) + 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 0, 1] = 1 / (Ca[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 1, 0] = 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 1, 1] = (r1 * e2 - r2 * e1) / math.sqrt(delta) - 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                B[k, 0] = 1 + ((r2 * e1 - r1 * e2) + 1 / (Ca[i] * Ra[i]) * (e1 - e2)) / math.sqrt(delta)
                B[k, 1] = 1 + (r2 * e1 - r1 * e2) / math.sqrt(delta)

        else:
            A = np.zeros((periods, 2, 2))
            B = np.zeros(periods, 2)
            r = -b / 2
            w = math.sqrt(-delta) / 2
            for k, t in enumerate(np.arange(0, 24, step)):
                A[k, 0, 0] = math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * (
                            1 / (Cm[i] * Rm[i]) + r)
                A[k, 0, 1] = math.exp(r * t) * math.sin(w * t) / w / (Ca[i] * Rm[i])
                A[k, 1, 0] = math.exp(r * t) * math.sin(w * t) / w / (Cm[i] * Rm[i])
                A[k, 1, 1] = math.exp(r * t) * math.cos(w * t) - math.exp(r * t) * math.sin(w * t) / w * (
                            1 / (Cm[i] * Rm[i]) + r)

                B[k, 0] = 1 - math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * (
                            1 / (Ca[i] * Ra[i]) + r)
                B[k, 1] = 1 - math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * r

        # 建立温度与功率之间的关系
        # Ta,Tm都是指t=[1,96],其余的是[0,95]
        for k in range(periods):
            if k == 0:
                Pac[i, k] = ((Ta[i, k] - A[k, 0, 0] * Ta0[i] - A[k, 0, 1] * Tm0[i]) / B[k, 0] - To[k])/(Cope[i] * Ra[i])
                if Pac[i, k] >= 0:
                    Pac[i, k] = 0
                Ta[i, k] = A[k, 0, 0] * Ta0[i] + A[k, 0, 1] * Tm0[i] + B[k, 0] * (To[k] + Pac[i, k] * Cope[i] * Ra[i])
                Tm[i, k] = A[k, 1, 0] * Ta0[i] + A[k, 1, 1] * Tm0[i] + B[k, 1] * (To[k] + Pac[i, k] * Cope[i] * Ra[i])
            else:
                Pac[i, k] = ((Ta[i, k] - A[k, 0, 0] * Ta[i, k - 1] - A[k, 0, 1] * Tm[i, k - 1]) / B[k, 0] - To[k]) / (Cope[i] * Ra[i])
                if Pac[i, k] >= 0:
                    Pac[i, k] = 0
                Ta[i, k] = A[k, 0, 0] * Ta[i, k - 1] + A[k, 0, 1] * Tm[i, k - 1] + B[k, 0] * (To[k] + Pac[i, k] * Cope[i] * Ra[i])
                Tm[i, k] = A[k, 1, 0] * Ta[i, k - 1] + A[k, 1, 1] * Tm[i, k - 1] + B[k, 1] * (To[k] + Pac[i, k] * Cope[i] * Ra[i])
            Pac[i, k] = -Pac[i, k]

    # for i in range(n_AC):
    #     plt.plot(range(periods),Ta[i])
    #     plt.plot(range(periods), Pac[i])
    # plt.plot(range(periods),To)
    # plt.show()

    for i in range(n_AC):
        ACParam[f'AC_{i}'] = {
            # 'To': To.tolist(),
            'Ta0': Ta0[i],
            'Tm0': Tm0[i],
            'Ca': Ca[i],
            'Cm': Cm[i],
            'Ra': Ra[i],
            'Rm': Rm[i],
            'Cope': Cope[i],
            'EnergyIni': sum(Pac[i]*24/periods)
        }


    return Pac, ACParam,To

# for i in range(n_AC):
#     plt.plot(range(periods),Ta[i])
#     plt.plot(range(periods), Pac[i])
# # plt.plot(range(periods),To)
# plt.show()
'''----------------------------------------------------------------------------------------------------------'''
n_EV = 50
SocMin = 0
SocMax = 80
StayMin = 10 #h
StayMax = 20 #h
Capicity = 80
ChargePowerMax = 7
Ev_n = 0.9
SOCDemand = 90
'''生成用户返回时间'''
def ComingTimePdf(x):
    if 0 <= x < 5.65:
        return (
                1/3.66/np.sqrt(2*np.pi)*
                np.exp(-(x+24-17.65)**2/2/3.66**2)
        )
    else:
        return (
                1/3.66/np.sqrt(2*np.pi)*
                np.exp(-(x-17.65)**2/2/3.66**2)
        )
# 计算概率密度函数对应的累积分布函数
def ComingTimeCdf(x):
    integral, _ = quad(ComingTimePdf, 0, x)
    return integral
# 计算 CDF 的逆函数
def ComingTime_inverse_cdf(p):
    x_values = np.linspace(0, 24, 1000)
    cdf_values = np.array([ComingTimeCdf(x) for x in x_values])
    return np.interp(p, cdf_values / cdf_values[-1], x_values)

'''生成用户停驻时长'''
def StayingTimePdf(x):
    if 0 <= x < 13.93-12:
        return (
                1/4.84/np.sqrt(2*np.pi)*
                np.exp(-(x+24-13.93)**2/2/4.84**2)
        )
    else:
        return (
                1/4.84/np.sqrt(2*np.pi)*
                np.exp(-(x-13.93)**2/2/4.84**2)
        )

def StayingTimeCdf(x):
    integral, _ = quad(StayingTimePdf, 0, x)
    return integral

    # 计算 CDF 的逆函数
def StayingTime_inverse_cdf(p, StayMin, StayMax):
    x_values = np.linspace(StayMin, StayMax, 1000)
    cdf_values = np.array([StayingTimeCdf(x) for x in x_values])
    return np.interp(p, cdf_values / cdf_values[-1], x_values)


'''生成用户初始SOC'''
def SOCIniPdf(x):
    if x == 0:
        return 0
    else:
        return (
                1 / 14.745 / np.sqrt(2 * np.pi) *
                np.exp(-(x - 57.64) ** 2 / 2 / 14.745 ** 2)
        )
def SOCIniCdf(x):
    integral, _ = quad(SOCIniPdf, 0, x)
    return integral
# 计算 CDF 的逆函数
def SOCIni_inverse_cdf(p, SocMin,SocMax):
    x_values = np.linspace(SocMin, SocMax, 1000)
    cdf_values = np.array([SOCIniCdf(x) for x in x_values])
    return np.interp(p, cdf_values / cdf_values[-1], x_values)

def EVDataGenerator(periods,n_EV,StayMin,StayMax,SocMin,SocMax,ChargePowerMax,SOCDemand,Ev_n,random_seed=40):
    EVsPower = np.zeros((n_EV, periods)) #存储原始功率
    EVsParam = {}                       #存储各个EV功率
    # 设置随机数种子以确保结果可复现
    np.random.seed(random_seed)
    # 生成均匀分布的随机数（0 到 1 之间的随机数）
    random_uniform = np.random.rand(n_EV)

    # 使用逆函数将均匀分布的随机数映射到目标分布上
    EVsComingTime_samples = ComingTime_inverse_cdf(random_uniform)
    EVsStayingTime_samples = StayingTime_inverse_cdf(random_uniform,StayMin,StayMax)
    EVsSocIni_samples = SOCIni_inverse_cdf(random_uniform,SocMin,SocMax)
    EVsDataList = {f"EV_{i}": {"EVsComingTime":int(EVsComingTime_samples[i]*periods/24),
                               "EVsStayingTime":int(EVsStayingTime_samples[i]*periods/24),
                               "EVsSocIni":EVsSocIni_samples[i],
                               }
                                 for i in range(n_EV)}

    # # 创建一个包含3个子图的图形，排成一行
    # fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    # # 可视化结果
    # axes[0].hist(EVsComingTime_samples, bins=30, density=True, alpha=0.7, label='Generated Samples')
    # x_values = np.linspace(0, 24, 1000)
    # pdf_values = np.array([ComingTimePdf(x) for x in x_values])
    # axes[0].plot(x_values, pdf_values, label='True PDF')
    # axes[0].set_xlabel('Values')
    # axes[0].set_ylabel('Density')
    # axes[0].set_title('Generated EVsComingTime Samples with True PDF')
    # axes[0].legend()
    #
    #
    # axes[1].hist(EVsStayingTime_samples, bins=30, density=True, alpha=0.7, label='Generated Samples')
    # x_values = np.linspace(0, 24, 1000)
    # pdf_values = np.array([StayingTimePdf(x) for x in x_values])
    # axes[1].plot(x_values, pdf_values, label='True PDF')
    # axes[1].set_xlabel('Values')
    # axes[1].set_ylabel('Density')
    # axes[1].set_title('Generated EVsStayingTime Samples with True PDF')
    # axes[1].legend()
    #
    #
    # axes[2].hist(EVsSocIni_samples, bins=30, density=True, alpha=0.7, label='Generated Samples')
    # x_values = np.linspace(0, SocMax, 1000)
    # pdf_values = np.array([SOCIniPdf(x) for x in x_values])
    # axes[2].plot(x_values, pdf_values, label='True PDF')
    # axes[2].set_xlabel('Values')
    # axes[2].set_ylabel('Density')
    # axes[2].set_title('Generated SocIni Samples with True PDF')
    # axes[2].legend(loc='upper left')
    # # plt.legend()
    # # 调整布局
    # plt.tight_layout()
    # plt.show()

    '''-------------------------------------开始生成充电功率序列----------------------------------'''
    for i in range(n_EV):
        EVsComingTime = EVsDataList[f'EV_{i}']['EVsComingTime']
        EVsStayingTime = EVsDataList[f'EV_{i}']['EVsStayingTime']
        EVsSocIni = EVsDataList[f'EV_{i}']['EVsSocIni']

        ChargingStartTime = EVsComingTime
        # print("ChargingStartTime:",ChargingStartTime)
        ChargingEndTime = (EVsComingTime + EVsStayingTime) % periods
        EnergyDemand = (SOCDemand - EVsSocIni)/100*Capicity   #电量Kwh
        ChargingTime = int(EnergyDemand/(ChargePowerMax*Ev_n*(24/periods)))+1 # +1保证充到EnergyDemand
        if ChargingTime > EVsStayingTime:
            print("该用户无法达到90%!")
            # break
            ChargingTime = EVsStayingTime
        # print("ChargingTime:",ChargingTime)
        if ChargingStartTime + ChargingTime <= periods:
            for T in range(ChargingStartTime,ChargingStartTime + ChargingTime):
                EVsPower[i,T] = ChargePowerMax
        else:
            for T in range(ChargingStartTime,periods):
                EVsPower[i,T] = ChargePowerMax
            for T in range(ChargingStartTime + ChargingTime - periods):
                EVsPower[i,T] = ChargePowerMax
        EVsParam[f'EV_{i}'] = {
                                "ChargingStartTime": ChargingStartTime,
                                "ChargingEndTime": ChargingEndTime,
                                "EVsSocIni": EVsSocIni,
                                'EnergyIni': sum(EVsPower[i]*Ev_n*24/periods),  # EV充电量
                                'ChargePowerMax':ChargePowerMax,
                                'Ev_n': Ev_n
                               }

    return EVsPower,EVsParam
    # plt.plot(range(len(IniPower)), IniPower,label = season[n])
    # plt.legend()
    # plt.show()
'''----------------------------------------------------
------------------------------------------------------------'''
def TFLPowerModel(periods,Pmin,Pmax):
    #生成随机负荷功率序列（假设范围在0到10之间，单位为千瓦）
    charging_power_sequence = np.round(np.random.uniform(Pmin, Pmax, periods),3)
    return charging_power_sequence.tolist()
'''生成一组可转移负荷数据'''
def TFLDataGenerator(periods,n_TFL,Pmin,Pmax,seed):
        TFLsPower = np.zeros((n_TFL,periods))
        TFLsParam = {}
        np.random.seed(seed)
        for i in range(n_TFL):
            TFLsPower[i] = TFLPowerModel(periods,Pmin,Pmax)
            TFLsParam['TFL_{i}'] = {'EnergyIni':sum(TFLsPower[i])*24/periods,
                                    'Pmax': Pmax}
        return TFLsPower,TFLsParam

def TimeShotConvert(data,periods):
    interval = round(96/periods)
    if interval == 1:
        return data
    else:
        NewData = data.reshape(-1, interval).mean(axis=1)
        return NewData

def PVDataGenerator(periods,PV_k):
    PVPower_json_file_path = 'PVPower.json'
    with open(PVPower_json_file_path, 'r') as json_file:
        PV = json.load(json_file)
    PVPower = TimeShotConvert(np.array(PV['summer']), periods) * PV_k
    return PVPower

def LoadDataGenerator(periods, LD_k):
    LDPower_json_file_path = 'LoadPower.json'
    with open(LDPower_json_file_path, 'r') as json_file:
        LD = json.load(json_file)
    LoadPower = TimeShotConvert(np.array(LD['summer']), periods) * LD_k
    return LoadPower

if __name__ == '__main__':
        periods = 24
        n_EV = 50
        SocMin = 0
        SocMax = 80
        StayMin = 10  # h
        StayMax = 20  # h
        Capicity = 80
        ChargePowerMax = 7
        Ev_n = 0.9
        SOCDemand = 90
        # '''生成TL所需参数'''
        # n_TFL = 200
        # P_max = 1
        # P_min = 5
        '''生成PV与LD所需参数'''
        PV_k = 1
        LD_k = 0.1

        AreaID = ['Area 1', 'Area 2', 'Area 3']
        Params = {} #存放各个台区的资源参数最终写入一个json文件中
        DERsPowerList = [] #存放各个台区的资源功率序列，并最终将其合并
        PowerRowIndexList = []
        '''生成AC调节温度上下限与调节模式'''
        Params['Mode'] = 1  # 调节模式：1 ：AC+EV 2：AC 3：EV
        Params['AC_Tmax'] = 28
        Params['AC_Tmin'] = 24
        '''依次生成各个台区的数据'''
        for s, area in enumerate(AreaID):
                n_AC = 300
                n_EV = 50
                seed = 40 + s
                ACPower,ACParam,To = ACPowerGenerator(periods,n_AC,seed)
                EVPower,EVParam = EVDataGenerator(periods,n_EV,StayMin,StayMax,SocMin,SocMax,ChargePowerMax,SOCDemand,Ev_n,seed)
                PVPower = PVDataGenerator(periods,PV_k)
                LDPower = LoadDataGenerator(periods, LD_k)

                # TFLPower,TFLParam = TFLDataGenerator(periods,n_TFL,P_min,P_max)
                print(area + '生成完毕！')
                Params[area] = {
                    'ACsParam': ACParam,
                    'EVParam': EVParam,
                    # 'PVParam': PVPower.tolist(),
                    # 'LDParam': LDPower.tolist()
                }
                Power = np.vstack((ACPower,EVPower))
                Power = np.vstack((Power,PVPower))
                Power = np.vstack((Power,LDPower))
                Power = np.vstack((Power,To))
                DERsPowerList.append(Power)
                PowerRowIndexList.append([area + '_' + f'AC_{i}' + '_Ini' for i in range(n_AC)]
                                     + [area + '_' + f'EV_{i}' + '_Ini' for i in range(n_EV)]
                                    + [area + '_' + 'PV'] + [area + '_' + 'LD'] + [area + '_' + 'Temperature'])

        DERsPower = np.vstack(DERsPowerList)  # 按顺序存放每一个资源的初始功率曲线
        PowerRowIndex = sum(PowerRowIndexList, [])  # 生成行名称，即每一个资源的名称
        # 生成时间列名称，从2024-10-18 00:00开始，间隔1小时
        start_time = datetime(2024, 10, 18, 0, 0)
        time_intervals = [start_time + timedelta(hours=i) for i in range(24)]
        TimeColumnIndex = [t.strftime('%Y-%m-%d %H:%M') for t in time_intervals]  # 列名称为时间格式

        # 创建DataFrame
        df = pd.DataFrame(DERsPower, index=PowerRowIndex, columns=TimeColumnIndex)
        df = df.T
        print("功率的DataFrame创建完毕！")
        csv_file_path = 'IniPower.csv'
        df.to_csv(csv_file_path, index=True)  # 保存到指定路径
        # tools.display_dataframe_to_user(name="带时间列名称的DataFrame", dataframe=df_time)

        json_file_path = 'DERParams.json'
        with open(json_file_path, 'w') as json_file:
            json.dump(Params, json_file, indent=4)
        print("参数已写入json中！")





